/*
 * CPLEXProblemSolver.c
 *
 *  Created on: 28/apr/2013
 *      Author: Ialenti Andrea, Marangon Fabio
 */

#include "../headers/Callback.h"
#include "../headers/CPLEXProblemGenerator.h"
#include "../headers/CPLEXProblemSolver.h"
#include "../headers/CPLEXUtility.h"
#include "../headers/Graph.h"
#include "../headers/Heuristic.h"
#include "../headers/Lagrangian.h"
#include "../headers/OneTree.h"
#include "../headers/Printer.h"
#include "../headers/StatusCodes.h"
#include "../headers/Utility.h"

#include <math.h>
#include <stdio.h>
#include <string.h>

int findSECs(const Config *cfg, CPXENVptr env, CPXLPptr lp, const Problem *p, SEC *secs, int *totSECs, bool cb, cbData *myData, void *cbdata, int from, bool verbose) {
	char *funcName = initCharArray(cfg->NAME_LEN);
	sprintf(funcName, "findSECs");

	int i, j, idx1, idx2, offset = 0;
	int nCycles = 0;
	int nEdges, cpxStat;
	Status status;
	double *x;
	int **adjMat, *cycLengths, *cycles;
	int nodes = cb ? myData->p->nodes : p->nodes;

	//	Allocate the required memory.
	if (verbose)
		fprintf(stdout, "\t\t\t'%s': allocating required memory. . .\n", funcName);
	adjMat = initAdjMat(nodes);
	cycLengths = initIntArray(nodes / 3);
	cycles = initIntArray(nodes);
	if (adjMat == NULL || cycLengths == NULL || cycles == NULL) {
		fprintf(stdout, "***'%s': memory allocation error. . .\n", funcName);
		freeCharArray(funcName);
		if (adjMat != NULL)
			freeAdjMat(adjMat, nodes);
		if (cycLengths != NULL)
			freeIntArray(cycLengths);
		if (cycles != NULL)
			freeIntArray(cycles);
		return getStatusCode(ERR_MEM_ALLOC);
	}
	if (cb) {
		nEdges = myData->numcols;
		x = initDblArray(nEdges);
		if (x == NULL) {
			fprintf(stdout, "***'%s': memory allocation error. . .\n", funcName);
			freeCharArray(funcName);
			freeAdjMat(adjMat, nodes);
			freeIntArray(cycLengths);
			freeIntArray(cycles);
			return getStatusCode(ERR_MEM_ALLOC);
		}
	}
	if (verbose)
		fprintf(stdout, "\t\t\t'%s': memory allocation completed. . .\n", funcName);

	//	Get the list of edges chosen in the current CPLEX solution.
	if (cb) {
		cpxStat = CPXgetcallbacknodex(env, cbdata, from, x, 0, nEdges - 1);
		if (cpxStat != CPX_SUCCESS) {
			fprintf(stdout, "***'%s': error retrieving the solution 'x'. CPLEX status: %d. . .\n", funcName, cpxStat);
			freeCharArray(funcName);
			freeAdjMat(adjMat, nodes);
			freeIntArray(cycLengths);
			freeIntArray(cycles);
			freeDblArray(x);
			return getStatusCode(ERR_CPX);
		}
	}
	else {
		status = getCPLEXSol(cfg, env, lp, &x, &nEdges);
		if (status != SUCCESS) {
			fprintf(stdout, "***'%s': error retrieving the solution 'x'. . .\n", funcName);
			freeCharArray(funcName);
			freeAdjMat(adjMat, nodes);
			freeIntArray(cycLengths);
			freeIntArray(cycles);
			return getStatusCode(ERR_GENERIC);
		}
	}

	//	Create the nodes adjacency matrix, linking each node with the two nodes that are connected in the current CPLEX solution.
	if (verbose)
		fprintf(stdout, "\t\t\t'%s': creating the nodes adjacency matrix. . .\n", funcName);
	if (cb)
		myData->totalCost = 0;
	resetAdjMatVals(adjMat, nodes);
	for (i = 0, j = 0; i < nEdges; i++)
		if (x[i] > 0.9) {

			if (cb)
				getEdgeNodes(myData->p, i, &idx1, &idx2);
			else
				getEdgeNodes(p, i, &idx1, &idx2);

			//	Save the 'idx1' adjacent.
			if (adjMat[idx1][0] == -1)
				adjMat[idx1][0] = idx2;
			else
				adjMat[idx1][1] = idx2;

			//	Save the 'idx2' adjacent.
			if (adjMat[idx2][0] == -1)
				adjMat[idx2][0] = idx1;
			else
				adjMat[idx2][1] = idx1;

			if (cb)
				myData->totalCost += myData->p->XeCost[i];
		}

	//	Find all the sub-tours in O(n²).
	if(verbose)
		fprintf(stdout, "\t\t\t'%s': searching all the sub-tours (O(n²)). . .\n", funcName);
	status = visitCycle(cfg, adjMat, nodes, cycLengths, &nCycles, cycles);
	if (status != SUCCESS) {
		fprintf(stdout, "***'%s': error searching for sub-tours. . .\n", funcName);
		freeCharArray(funcName);
		freeAdjMat(adjMat, nodes);
		freeIntArray(cycLengths);
		freeIntArray(cycles);
		freeDblArray(x);
		return getStatusCode(ERR_GENERIC);
	}

	//	Create one SEC for every tour found.
	if (verbose)
		fprintf(stdout, "\n\t\t\t'%s': creating SECs for the sub-tours found. . .\n\n", funcName);
	if (nCycles != 1)
		for (i = 0, offset = 0; i < nCycles; i++, offset += cycLengths[i - 1]) {

			//	Set the right hand side for the current SEC, and the number of non zero coefficients.
			secs[i].rhs[0] = cycLengths[i] - 1;
			secs[i].sense[0] = 'L';
			secs[i].matbeg[0] = 0;
			secs[i].nonZero = cycLengths[i];

			//	Set the new constraint name.
			sprintf(secs[i].constrName[0], "c_SEC_%d", (*totSECs)++);

			//	Set the left hand side of the constraint.
			for (j = 0; j < secs[i].nonZero - 1; j++) {
				if (cb)
					secs[i].matind[j] = getEdgeIndex(myData->p, cycles[offset + j], cycles[offset + j + 1]);
				else
					secs[i].matind[j] = getEdgeIndex(p, cycles[offset + j], cycles[offset + j + 1]);
				secs[i].matval[j] = 1.0;
			}
			if (cb)
				secs[i].matind[secs[i].nonZero - 1] = getEdgeIndex(myData->p, cycles[offset + secs[i].nonZero - 1], cycles[offset]);
			else
				secs[i].matind[secs[i].nonZero - 1] = getEdgeIndex(p, cycles[offset + secs[i].nonZero - 1], cycles[offset]);
			secs[i].matval[secs[i].nonZero - 1] = 1.0;
		}

	//	Free memory and return the number of cycles found.
	freeCharArray(funcName);
	freeIntArray(cycLengths);
	freeIntArray(cycles);
	freeAdjMat(adjMat, nodes);
	freeDblArray(x);
	return nCycles;
}

Status visitCycle(const Config *cfg, int **adjMat, int nodes, int *cLengths, int *nCycles, int *cycles) {
	char *funcName = initCharArray(cfg->NAME_LEN);
	sprintf(funcName, "visitCycle");

	int i, curNode = 0, pos = 0, curCycLength = 0;
	bool unvisited = true;

	int *visited = initIntArray(nodes);
	if (visited == NULL) {
		fprintf(stdout, "***'%s': memory allocation error!\n", funcName);
		freeCharArray(funcName);
		return ERR_MEM_ALLOC;
	}
	*nCycles = 0;
	cycles[0] = curNode;

	//	Initialize all nodes to unvisited.
	for (i = 0; i < nodes; i++)
		visited[i] = 0;

	//	Repeat until all nodes have been visited.
	while (unvisited) {

		//	Flag the current node as visited.
		visited[curNode] = 1;

		//	Check if any of the node linked to the current one is still unvisited, if it is use it as next node.
		if (visited[adjMat[curNode][0]] == 0) {
			curNode = adjMat[curNode][0];
			curCycLength++;
			cycles[++pos] = curNode;
		} else if (visited[adjMat[curNode][1]] == 0) {
			curNode = adjMat[curNode][1];
			curCycLength++;
			cycles[++pos] = curNode;
		}

		//	If there are no unvisited nodes linked to the current one than the cycle is complete.
		else {

			//	Save the current cycle length and update the cycles counter.
			cLengths[*nCycles] = curCycLength + 1;
			(*nCycles)++;

			//	Search for a new unvisited node to start from to look for another cycle.
			for (i = 0, unvisited = false, curCycLength = 0; i < nodes; i++)
				if (visited[i] == 0) {
					curNode = i;
					cycles[++pos] = curNode;
					unvisited = true;
					break;
				}
		}
	}

	//	Free memory.
	freeCharArray(funcName);
	freeIntArray(visited);
	return SUCCESS;
}

Status solveIterMiliotis(const Config *cfg, CPXENVptr env, CPXLPptr lp, const Problem *p, const Graph *g, bool verbose, bool ps, time_t t1) {
	char *funcName = initCharArray(cfg->NAME_LEN);
	char *probName = initCharArray(cfg->NAME_LEN);
	int i, cpxStat, solCnt = 1, totSECs = 0;
	Status status;
	sprintf(funcName, "solveIterMiliotis");
	sprintf(probName, "%s%s.lp", cfg->CPX_FILE_PATH, p->name);
	time_t t2;

	//	Set the global time limit depending on the problem size.
	int globalTimeLimit = cfg->CPX_PS_GTL;
	if (p->nodes > 400)
		globalTimeLimit *= 6;
	else if (p->nodes > 200)
		globalTimeLimit *= 4;
	else if (p->nodes > 100)
		globalTimeLimit *= 2;

	//	Allocate an array of SECs that can store the maximum number of sub-tour for the problem.
	SEC *secs;
	if (!initSECsArray(cfg, &secs, p->nodes / 3, p->nodes)) {
		fprintf(stdout, "***'%s': SECs memory allocation error. . .\n", funcName);
		freeCharArray(funcName);
		freeCharArray(probName);
		return ERR_MEM_ALLOC;
	}

	//	Check if the first solution is an Hamiltonian cycle (total number of cycles is 1).
	if (verbose)
		fprintf(stdout, "\t\t'%s': searching for tours in the first solution. . .\n", funcName);
	int cycles = findSECs(cfg, env, lp, p, secs, &totSECs, false, NULL, NULL, 0, verbose);

	//	If any SEC was found, add them to the problem.
	if (cycles > 1) {
		for (i = 0; i < cycles; i++) {
			cpxStat = CPXaddrows(env, lp, secs[i].nonZero, 1, secs[i].nonZero, secs[i].rhs, secs[i].sense, secs[i].matbeg, secs[i].matind, secs[i].matval, NULL, secs[i].constrName);
			if (cpxStat != CPX_SUCCESS) {
				fprintf(stdout, "***'%s': error adding %d SECs to the problem using CPXaddrows. CPLEX status: %d. . .\n", funcName, cycles, cpxStat);
				freeCharArray(funcName);
				freeCharArray(probName);
				disposeSECsArray(secs, p->nodes / 3);
				return ERR_CPX;
			}
		}
		if (verbose)
			fprintf(stdout, "\t\t'%s': %d SECs added. . .\n", funcName, cycles);
	}
	else if (cycles < 0) {
		if (verbose)
			fprintf(stdout, "\t\t'%s': unable to identify any tour. . .\n", funcName);
		freeCharArray(funcName);
		freeCharArray(probName);
		disposeSECsArray(secs, p->nodes / 3);
		return ERR_GENERIC;
	}

	//	Iteratively add SECs and re-resolve the problem until a single Hamiltonian cycle is found.
	if (verbose)
		fprintf(stdout, "\n\t\t'%s': running iterative algorithm. . .\n", funcName);
	while (cycles != 1) {
		solCnt++;

		if (verbose)
			fprintf(stdout, "\t\t'%s': re-running MIP opt (#%d) with the newly added SECs. . .\n\n", funcName, solCnt);

		//	Check the time constraint.
		time(&t2);
		if(ps) {
			if(difftime(t2, t1) * cfg->NORM_TIME > globalTimeLimit)
				return ERR_TIME;
		}
		else {
			if(difftime(t2, t1) * cfg->NORM_TIME > cfg->CPX_TIME_LIM)
				return ERR_TIME;
		}

		//	Re-optimize the problem with the new SECs and plot the solution.
		cpxStat = CPXmipopt(env, lp);
		if (cpxStat != CPX_SUCCESS) {
			fprintf(stdout, "***'%s': failed to optimize MIP # %d. Status: %d.\n", funcName, solCnt, cpxStat);
			freeCharArray(funcName);
			freeCharArray(probName);
			disposeSECsArray(secs, p->nodes / 3);
			return ERR_CPX;
		}

		//	Plot the current solution.
		if (verbose) {
			double *x;
			int nEdges;
			status = getCPLEXSol(cfg, env, lp, &x, &nEdges);
			if (status == SUCCESS) {
				fprintf(stdout, "\n\t\t'%s': plotting MIP opt #%d solution. . .\n", funcName, solCnt);
				if (ps)
					status = plotSolution(cfg, env, lp, p, g, solCnt, ps, false);
				else if (solCnt % 15 == 0)
					status = plotSolution(cfg, env, lp, p, g, solCnt, false, false);
				if (status != SUCCESS)
					fprintf(stdout, "\t***'%s': error plotting CPLEX solution #%d. . .\n", funcName, solCnt);
			}
		}

		//	Search for SECs on the new solution.
		if (verbose)
			fprintf(stdout, "\t\t'%s': searching for tours in the solution #%d. . .\n", funcName, solCnt);
		cycles = findSECs(cfg, env, lp, p, secs, &totSECs, false, NULL, NULL, 0, verbose);
		if (cycles > 1) {
			for (i = 0; i < cycles; i++) {
				cpxStat = CPXaddrows(env, lp, secs[i].nonZero, 1, secs[i].nonZero, secs[i].rhs, secs[i].sense, secs[i].matbeg, secs[i].matind, secs[i].matval, NULL, secs[i].constrName);
				if (cpxStat != CPX_SUCCESS) {
					fprintf(stdout, "***'%s': error adding %d SECs to the problem using CPXaddrows. CPLEX status: %d. . .\n", funcName, cycles, cpxStat);
					freeCharArray(funcName);
					freeCharArray(probName);
					disposeSECsArray(secs, p->nodes / 3);
					return ERR_CPX;
				}
			}
			if (verbose)
				fprintf(stdout, "\t\t'%s': %d SECs added. . .\n", funcName, cycles);
		}
		else if (cycles < 0) {
			if (verbose)
				fprintf(stdout, "***'%s': unable to identify any tour. . .\n", funcName);
			freeCharArray(funcName);
			freeCharArray(probName);
			disposeSECsArray(secs, p->nodes / 3);
			return ERR_GENERIC;
		}
	}
	if (verbose)
		fprintf(stdout, "\n\t\t'%s': iterative algorithm completed. . .\n\n", funcName);

	//	Write the final problem with all the necessary SECs to file.
	cpxStat = CPXwriteprob(env, lp, probName, "lp");
	if (cpxStat != CPX_SUCCESS) {
		fprintf(stdout, "\n***'%s': error writing problem to file. CPLEX status: %d.\n", funcName, cpxStat);
		freeCharArray(funcName);
		freeCharArray(probName);
		disposeSECsArray(secs, p->nodes / 3);
		return ERR_CPX;
	}

	//	Free memory.
	freeCharArray(funcName);
	disposeSECsArray(secs, p->nodes / 3);
	freeCharArray(probName);
	return SUCCESS;
}

Status proxySearch(const Config *cfg, CPXENVptr env, CPXLPptr lp, CPXLPptr psLP, Problem *p, const Problem *psP, Graph *g, bool verbose) {
	Status status;
	time_t t1, t2;
	cbData data;
	int cpxStat, i, idx1, idx2, nEdges, cnt = 0;
	double *x, solCost;
	char *warmStart = initCharArray(cfg->NAME_LEN);
	char *title = initCharArray(cfg->NAME_LEN);
	char *file = initCharArray(cfg->NAME_LEN);
	char *funcName = initCharArray(cfg->NAME_LEN);
	sprintf(funcName, "proxySearch");
	sprintf(warmStart, "%swarmStart_%s.mst", cfg->CPX_FILE_PATH, psP->name);

	//	Set the global time limit depending on the problem size.
	int globalTimeLimit = cfg->CPX_PS_GTL;
	if (p->nodes > 400)
		globalTimeLimit *= 6;
	else if (p->nodes > 200)
		globalTimeLimit *= 4;
	else if (p->nodes > 100)
		globalTimeLimit *= 2;

	//	Set CPLEX to stop when the first integer solution is found.
	cpxStat = CPXsetintparam(env, CPX_PARAM_INTSOLLIM, 1);
	if (cpxStat != CPX_SUCCESS) {
		fprintf(stdout, "***'%s': failed to set 'CPX_PARAM_INTSOLLIM'.\n", funcName);
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}

	//	Set tree memory size limit and execution time limit for CPLEX.
	cpxStat = CPXsetdblparam(env, CPX_PARAM_TRELIM, cfg->CPX_PS_ML);
	if (cpxStat != CPX_SUCCESS) {
		fprintf(stdout, "***'%s': failed to set 'CPX_PARAM_TRELIM'.\n", funcName);
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}
	cpxStat = CPXsetdblparam(env, CPX_PARAM_TILIM, cfg->CPX_PS_TL / cfg->NORM_TIME);
	if (cpxStat != CPX_SUCCESS) {
		fprintf(stdout, "***'%s': failed to set 'CPX_PARAM_TILIM'.\n", funcName);
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}

	//	Find the first proximity search solution.
	if (verbose)
		fprintf(stdout, "\t'%s': running first proximity search MIP opt. . .\n", funcName);
	cpxStat = CPXmipopt(env, psLP);
	if (cpxStat != 0) {
		fprintf(stdout, "***'%s': failed to optimize MIP. Status: %d.\n", funcName, cpxStat);
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}
	int s = CPXgetstat(env, psLP);
	if (verbose)
		fprintf(stdout, "\t'%s': first MIP opt completed. . .\n", funcName);

	//	Start the proximity search iterative algorithm.
	time(&t1);
	do {
		cnt++;

		//	Remove all sub-tours from the current solution using Miliotis' method.
		if (verbose)
			fprintf(stdout, "\t\t'%s': MIP OPT status: %d. Removing sub-tours using Miliotis' method for MIP opt run #%d. . .\n", funcName, s, cnt);
		status = solveIterMiliotis(cfg, env, psLP, psP, g, false, true, t1);
		if (status == SUCCESS) {
			status = getCPLEXSol(cfg, env, psLP, &x, &nEdges);
			if (status != SUCCESS) {
				fprintf(stdout, "***'%s': failed to retrieve proximity search solution #%d.\n", funcName, cnt);
				freeCharArray(funcName);
				freeCharArray(warmStart);
				freeCharArray(title);
				freeCharArray(file);
				return ERR_GENERIC;
			}

			if(verbose) {
				time(&t2);
				for (i = 0, solCost = 0; i < psP->edges; i++)
					if (x[i] > 1 - cfg->EPS) {
						getEdgeNodes(psP, i, &idx1, &idx2);
						solCost += g->costs_origin[idx1][idx2];
					}
				fprintf(stdout, "\t\t'%s': sub-tours removed. . .\n", funcName);
				fprintf(stdout, "\t\t'%s': current solution cost: %0.2f. . .\n", funcName, solCost);
				fprintf(stdout, "\t\t'%s': elapsed time: %0.2f second(s) (%0.2f normalized time second(s)). . .\n", funcName, difftime(t2, t1), difftime(t2, t1) * cfg->NORM_TIME);
				fflush(stdout);
			}

			//	Plot 1 proximity search solution out of 3.
			if (cnt % 3 == 0) {
				status = plotSolution(cfg, env, psLP, psP, g, cnt, true, false);
				if (status != SUCCESS)
					fprintf(stdout, "\n***'%s': error plotting proximity search solution #%d. . .\n\n", funcName, cnt);
			}

			//	Store the last proximity search solution to be used as the warm start for the normal TSP problem.
			cpxStat = CPXwritemipstarts(env, psLP, warmStart, 0, 0);
			if (cpxStat != CPX_SUCCESS) {
				fprintf(stdout, "***'%s': failed to store proximity search solution #%d for warm start.\n", funcName, cnt);
				freeCharArray(funcName);
				freeCharArray(warmStart);
				freeCharArray(title);
				freeCharArray(file);
				return ERR_CPX;
			}
		}
		else if(status == ERR_TIME) {
			if (verbose) {
				time(&t2);
				fprintf(stdout, "\t\t'%s': global time limit reached. . . \n", funcName);
				fprintf(stdout, "\t\t'%s': elapsed time: %0.2f second(s) (%0.2f normalized time second(s)). . .\n", funcName, difftime(t2, t1), difftime(t2, t1) * cfg->NORM_TIME);
			}
			break;
		}
		else {
			if (verbose) {
				time(&t2);
				fprintf(stdout, "\t\t'%s': sub-tours removal error. . .\n", funcName);
				fprintf(stdout, "\t\t'%s': elapsed time: %0.2f second(s) (%0.2f normalized time second(s)). . .\n", funcName, difftime(t2, t1), difftime(t2, t1) * cfg->NORM_TIME);
			}
			break;
		}

		//	Get the current solution x.
		status = getCPLEXSol(cfg, env, psLP, &x, &nEdges);
		if (status != SUCCESS) {
			if (verbose)
				fprintf(stdout, "\t'%s': error retrieving MIP opt solution #%d. . .\n", funcName, cnt);
			break;
		}

		//	Calculate new RHS and change proximity search constraint and objective function.
		status = changePSconstr(cfg, env, psLP, psP, g, x, 1);
		if (status != SUCCESS) {
			if (verbose)
				fprintf(stdout, "\t'%s': error changing proximity search constraint and objective function for MIP opt #%d. . .\n", funcName, cnt + 1);
			break;
		}

		//	Stop the algorithm if the global time limit is reached.
		time(&t2);
		if (difftime(t2, t1) * cfg->NORM_TIME > globalTimeLimit) {
			if (verbose)
				fprintf(stdout, "\n\t'%s': proximity search algorithm reached global time limit. . .\n", funcName);
			break;
		}

		//	Solve the new problem.
		if (verbose)
			fprintf(stdout, "\n\t'%s': running proximity search MIP opt #%d. . .\n", funcName, cnt + 1);
		cpxStat = CPXmipopt(env, psLP);
		s = CPXgetstat(env, psLP);
	} while (cpxStat == CPX_SUCCESS);
	if(verbose)
		fprintf(stdout, "\n\t'%s': proximity search algorithm ended. . .\n", funcName);

	//	Calculate the last solution cost.
	for (i = 0, solCost = 0; i < psP->edges; i++)
		if (x[i] > 1 - cfg->EPS) {
			getEdgeNodes(psP, i, &idx1, &idx2);
			solCost += g->costs_origin[idx1][idx2];
		}

	//	Plot the last solution.
	sprintf(title, "Proximity search final solution for '%s'. Solution cost = %0.2f", psP->name, solCost);
	sprintf(file, "%sPS_solution.dat", cfg->PLT_PATH);
	status = plotCPLEXSol(cfg, x, psP, g, file, title, cfg->COLOR_BLUE);
	if (status != SUCCESS)
		fprintf(stdout, "\n***'%s': error plotting proximity search final solution. . .\n\n", funcName);

	//	Set CPLEX parameters to optimize the normal TSP problem.
	cpxStat = CPXsetintparam(env, CPX_PARAM_INTSOLLIM, 2100000000);
	if (cpxStat != CPX_SUCCESS) {
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}
	cpxStat = CPXsetdblparam(env, CPX_PARAM_TRELIM, cfg->CPX_MEM_LIM);
	if (cpxStat != CPX_SUCCESS) {
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}
	cpxStat = CPXsetdblparam(env, CPX_PARAM_TILIM, cfg->CPX_TIME_LIM / cfg->NORM_TIME);
	if (cpxStat != CPX_SUCCESS) {
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}

	//	Configure CPLEX to enable warm start.
	cpxStat = CPXsetintparam(env, CPX_PARAM_ADVIND, 1);
	if (cpxStat != CPX_SUCCESS) {
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}

	//	Initialize callback data.
	initCBData(cfg, &data, g, p, &lp, verbose);
	data.bestLSSol = solCost;

	//	Set the lazy constraint callback.
	cpxStat = CPXsetlazyconstraintcallbackfunc(env, cbLazySEC, &data);
	if (cpxStat != CPX_SUCCESS) {
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}

	//	Set the min cut callback.
	status = CPXsetusercutcallbackfunc(env, cbMinCut, &data);
	if (cpxStat != CPX_SUCCESS) {
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}

	//	Set MIP opt warm start.
	cpxStat = CPXreadcopymipstarts(env, lp, warmStart);
	if (cpxStat != CPX_SUCCESS) {
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}

	//	Turn on CPLEX output to the standard output (console).
	if(verbose)
		fprintf(stdout, "\t'%s': turning on CPLEX screen indicator. . .\n", funcName);
	cpxStat = CPXsetintparam(env, CPX_PARAM_SCRIND, CPX_ON);
	if (cpxStat != CPX_SUCCESS)
		fprintf(stdout, "***'%s': failed to turn on screen indicator, error #: %d.\n", funcName, cpxStat);

	//	Solve the TSP problem using callbacks.
	if(verbose)
		fprintf(stdout, "\t'%s': solving original TSP problem with CPLEX warm start with solution cost: %0.2f. . .\n\n\n", funcName, solCost);
	cpxStat = CPXmipopt(env, lp);
	if (cpxStat != CPX_SUCCESS) {
		freeCharArray(funcName);
		freeCharArray(warmStart);
		freeCharArray(title);
		freeCharArray(file);
		return ERR_CPX;
	}

	//	Check the CPLEX mip opt status.
	bool opt = true;
	cpxStat = CPXgetstat(env, lp);
	char code[CPXMESSAGEBUFSIZE];
	CPXgetstatstring(env, cpxStat, code);
	fprintf(stdout, "\n\t'%s': optimization completed with status %d ('%s'). . .\n", funcName, cpxStat, code);
	if (cpxStat != CPXMIP_OPTIMAL && cpxStat != CPXMIP_OPTIMAL_TOL) {
		opt = false;
		if (verbose)
			fprintf(stdout, "\n\t'%s': *WARNING* MIP opt was stopped before it could certify the tour's optimality. MIP OPT status: %d. . .\n\n", funcName, cpxStat);
	}

	//	Plot the optimal TSP solution.
	status = plotSolution(cfg, env, lp, p, g, 0, true, opt);
	if (status != SUCCESS)
		fprintf(stdout, "\n***'%s': error plotting TSP optimal solution. . .\n\n", funcName);

	//	Free the memory.
	freeCharArray(funcName);
	freeCharArray(warmStart);
	freeCharArray(title);
	freeCharArray(file);
	return SUCCESS;
}

Status changePSconstr(const Config *cfg, CPXENVptr env, CPXLPptr psLP, const Problem *p, const Graph *g, double *x, double theta) {
	char *funcName = initCharArray(cfg->NAME_LEN);
	sprintf(funcName, "changePSconstr");
	int j, i, idx1, idx2, sz, cnt, cpxStat;

	//	Calculate the new right hand side.
	sz = p->edges;
	cnt = p->numcols;
	double newRHS = -theta;
	for (j = 0; j < sz; j++)
		if (x[j] > 0.9) {
			getEdgeNodes(p, j, &idx1, &idx2);
			newRHS += g->costs_origin[idx1][idx2];
		}
	if (newRHS == -theta) {
		fprintf(stdout, "***'%s': error calculating the new right hand side. . .\n", funcName);
		freeCharArray(funcName);
		return ERR_GENERIC;
	}

	//	Get Proximity search constraint index.
	cpxStat = CPXgetrowindex(env, psLP, cfg->CPX_PS_C_NAME, &i);
	if (cpxStat != CPX_SUCCESS) {
		fprintf(stdout, "***'%s': error retrieving proximity search constraint index. . .\n", funcName);
		freeCharArray(funcName);
		return ERR_CPX;
	}

	//	Change proximity search constraint's right hand side.
	cpxStat = CPXchgrhs(env, psLP, 1, &i, &newRHS);
	if (cpxStat != CPX_SUCCESS) {
		fprintf(stdout, "***'%s': error setting new right hand side for the proximity search constraint. . .\n", funcName);
		freeCharArray(funcName);
		return ERR_CPX;
	}

	//	Change objective function.
	int indices[cnt];
	double vals[cnt];
	for (j = 0; j < cnt; j++) {
		if (x[j] == 1)
			vals[j] = -1;
		if (x[j] == 0)
			vals[j] = 1;

		indices[j] = j;
	}
	cpxStat = CPXchgobj(env, psLP, cnt, indices, vals);
	if (cpxStat != CPX_SUCCESS) {
		fprintf(stdout, "***'%s': error updating the objective function. . .\n", funcName);
		freeCharArray(funcName);
		return ERR_CPX;
	}

	//	Free memory.
	freeCharArray(funcName);
	return SUCCESS;
}

Status addLagrangianSECs(const Config *cfg, CPXENVptr env, CPXLPptr lp, const Problem *p, OneTree *lagrTree) {
	char *funcName = initCharArray(cfg->NAME_LEN);
	sprintf(funcName, "addLagrangianSECs");

	int i, howManySECs = 0;
	SEC *secs;
	Status status;

	//	Memory allocation.
	int *actualTree = initIntArray(lagrTree->size);
	int *alreadyRoot = initIntArray(lagrTree->size);
	if (!initSECsArray(cfg, &secs, p->nodes, p->nodes)) {
		fprintf(stdout, "***'%s': SECs memory allocation error. . .\n", funcName);
		freeCharArray(funcName);
		freeIntArray(actualTree);
		freeIntArray(alreadyRoot);
		return ERR_MEM_ALLOC;
	}

	//	Array initialization.
	for (i = 0; i < lagrTree->size; i++)
		actualTree[i] = 1;

	//	Extract SECs and add them to the problem.
	status = extractSECsFromOneTree(cfg, &env, &lp, lagrTree, secs, p, actualTree, alreadyRoot, &howManySECs);
	if (status != SUCCESS) {
		fprintf(stdout, "***'%s': error extracting SECs. . .\n", funcName);
		freeCharArray(funcName);
		freeIntArray(actualTree);
		freeIntArray(alreadyRoot);
		disposeSEC(secs);
		return ERR_GENERIC;
	}

	//	Free memory.
	freeCharArray(funcName);
	freeIntArray(actualTree);
	freeIntArray(alreadyRoot);
	disposeSEC(secs);
	return SUCCESS;
}

Status extractSECsFromOneTree(const Config *cfg, CPXENVptr *env, CPXLPptr *lp, const OneTree *t, SEC *secs, const Problem *p, int *actualTree, int *alreadyRoot, int *howManySECs) {
	char *funcName = initCharArray(cfg->NAME_LEN);
	sprintf(funcName, "extractSECsFromOneTree");

	int i, nodeA, nodeB, tmpNodeA, tmpNodeB, aCnt, bCnt, max = -cfg->INF;
	Status status;

	//	Find the maximum cost edge.
	for (i = 0; i < t->size; i++)

		//	Check if the edge exists.
		if (max < t->edges[i].cost && t->edges[i].cost != 0) {
			tmpNodeA = t->edges[i].start->id;
			tmpNodeB = t->edges[i].end->id;

			//	Check if the node belongs to the 'actualTree'.
			if (actualTree[tmpNodeA] == 1 && actualTree[tmpNodeB] == 1) {
				max = t->edges[i].cost;
				nodeA = t->edges[i].start->id;
				nodeB = t->edges[i].end->id;
			}
		}

	//	Clear 'actualTree' and reset the counters
	for (i = aCnt = bCnt = 0; i < t->size; actualTree[i++] = 0);

	//	Find the subtree of 'nodeA' and extract its SECs.
	visitSubtree(t, nodeA, &aCnt, nodeB, actualTree, alreadyRoot);
	status = buildSEC(cfg, env, lp, actualTree, p, secs, howManySECs);
	if (status != SUCCESS) {
		fprintf(stdout, "***'%s': error extracting the SECs from sub-tree of node %d. . .\n", funcName, nodeA);
		freeCharArray(funcName);
		return ERR_GENERIC;
	}
	//else
		//*howManySECs = *howManySECs + 1;
	alreadyRoot[nodeA] = 1;

	//	Run the function recursively.
	if (aCnt > 3)
		extractSECsFromOneTree(cfg, env, lp, t, secs, p, actualTree, alreadyRoot, howManySECs);

	//	Clear actualTree
	for (i = 0; i < t->size; actualTree[i++] = 0);

	//	Find the subtree of 'nodeB' and extract its SECs.
	visitSubtree(t, nodeB, &bCnt, nodeA, actualTree, alreadyRoot);
	status = buildSEC(cfg, env, lp, actualTree, p, secs, howManySECs);
	if (status != SUCCESS) {
		fprintf(stdout, "***'%s': error extracting the SECs from sub-tree of node %d. . .\n", funcName, nodeB);
		freeCharArray(funcName);
		return ERR_GENERIC;
	}
	//else
		//*howManySECs = *howManySECs + 1;
	alreadyRoot[nodeB] = 1;

	//	Run the function recursively.
	if (bCnt > 3)
		extractSECsFromOneTree(cfg, env, lp, t, secs, p, actualTree, alreadyRoot, howManySECs);
	return SUCCESS;
}

Status buildSEC(const Config *cfg, CPXENVptr *env, CPXLPptr *lp, int *nodes, const Problem *p, SEC *secs, int *howManySECs) {
	char *funcName = initCharArray(cfg->NAME_LEN);
	sprintf(funcName, "extractSECFromNodes");

	int edge, nonZero, matbeg[1], i, j, q;
	double rhs[1];
	char sense[1];

	int *edges = initIntArray(p->edges);
	if(edges == NULL) {
		fprintf(stdout, "***'%s': memory allocation error. . .\n", funcName);
		freeCharArray(funcName);
		return ERR_MEM_ALLOC;
	}
	double *vals = initDblArray(p->edges);
	if(vals == NULL) {
		fprintf(stdout, "***'%s': memory allocation error. . .\n", funcName);
		freeCharArray(funcName);
		freeIntArray(edges);
		return ERR_MEM_ALLOC;
	}

	//	Store the graph edges in the 'edges' array.
	for (i = 0, q = 0; i < p->nodes - 1; i++)
		if (nodes[i] == 1)
			for (j = i + 1; j < p->nodes; j++)
				if (nodes[j] == 1)
					edge = getEdgeIndex(p, i, j);
					if (edge >= 0)
						edges[q++] = edge;

	//	The minimum number of edges to build a SEC is 3.
	if (q >= 3) {

		//	Build the SEC.
		nonZero = q;
		rhs[0] = q - 1;
		sense[0] = 'L';
		matbeg[0] = 0;

		for (i = 0; i < q; i++)
			vals[i] = 1.0;

		//	Add the SEC to the problem.
		int cpxStat = CPXaddrows(*env, *lp, nonZero, 1, nonZero, rhs, sense, matbeg, edges, vals, NULL, NULL);
		if (cpxStat != CPX_SUCCESS) {
			fprintf(stdout, "***'%s': error adding the SEC. CPLEX status: %d.\n", funcName, cpxStat);
			freeCharArray(funcName);
			freeIntArray(edges);
			freeDblArray(vals);
			return ERR_CPX;
		}
		(*howManySECs)++;

		//	Free memory.
		freeCharArray(funcName);
		freeIntArray(edges);
		freeDblArray(vals);
		return SUCCESS;
	}

	//	Free memory.
	freeCharArray(funcName);
	freeIntArray(edges);
	freeDblArray(vals);
	return SUCCESS;
}

void visitSubtree(const OneTree *t, int root, int *actualNodeNumber, int father, int *actualTree, int *alreadyRoot) {
	int i;

	if (alreadyRoot[root] != 1) {
		actualTree[root] = 1;
		*actualNodeNumber = *actualNodeNumber + 1;

		for (i = 1; i < t->size; i++)
			if (t->costs[root][i] != 0 && i != father)
				visitSubtree(t, i, actualNodeNumber, root, actualTree, alreadyRoot);
	}
	return;
}
